Charged particles on a 2D plane subject to anisotropic Jahn- Teller interactions. 
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The properties of a system of charged particles on a 2D lattice, subject to an anisotropic Jahn- 
Teller-type interaction and 3D Coulomb repulsion are investigated. In the mean-field approximation 
without Coulomb interaction, the system displays a phase transition of first order. When the long 
range Coulomb interaction is included, Monte Carlo simulations show that the system displays very 
diverse mesoscopic textures, ranging from spatially disordered pairs to ordered arrays of stripes, 
or charged clusters, depending only on the ratio of the two interactions (and the particle density). 
Remarkably, charged objects with even number of particles are more stable than with odd number of 
particles. We suggest that the diverse functional behaviour - including superconductivity - observed 
in oxides can be thought to arise from the self-organization of this type. 
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The standard theoretical models of strongly corre- 
lated electrons, such as the Hubbard model[l| or the 
t — J modelQ neglect two important interactions, namely 
long-range Coulomb repulsion and lattice distortions 
caused by charged particles. Moreover, these quantum 
mechanical models are typically used to study T w 
properties. As such, these models have found limited 
applicability in predicting the finite-temperature func- 
tional behaviour in systems such as cuprate superconduc- 
tors and other oxides. An important aspect of the prob- 
lem which has been of great interest recently is the exis- 
tence of intrinsic mesoscale inhomogeneity in these sys- 
tems, for which there is mounting experimental evidence 
from neutron scattering^, XAFS 3 , STM|i and time- 
resolved carrier dynamics [6| amongst others |2|. There is 
emerging consensus that in doped cuprates charge car- 
riers may phase segregate to form nano-scale textures. 
These are believed to be of importance for achieving their 
functional properties, and particularly superconductiv- 
ity. The idea of charge segregation in cuprates appeared 
soon after the discovery of superconductivity [g, 0, El i 
but in most cases, long-range Coulomb repulsion was not 
considered. More recently it was suggested that interplay 
of short range lattice attraction and long-range Coulomb 
repulsion could lead to the formation short metallic or 
insulating strings of polarons[llL H^|. Since an isotropic 
interaction cannot lead to stripe formation we suggested 
instead that an anisotropic mesoscopic Jahn- Teller in- 
teraction between electrons and k ^ optical phonons 
might lead to the formation of pairs and stripesjlij. A 
slightly different approach, based on elasticity was con- 
sidered more recently for the case of manganites by Kugel 
and Khomskii u sing the methods of Eremin et al.|lj|, 
and by Shenoy et al. |lfj . The importance of the interplay 
of long-range and short range forces within an Ising-like 
model was discussed by Low et al.[r^. 

The fundamental question which we try and answer 
here is how charged particles order in the presence 
of anisotropic Jahn- Teller type interaction, particularly 
when their density becomes large. We consider charged 



particles on a 2D square lattice subject to only the 
long-range Coulomb interaction and an anisotropic Jahn- 
Teller (JT) deformation. In the mean field (MF) ap- 
proximation without Coulomb repulsion, the system dis- 
plays a first order phase transition to an ordered state 
below some critical temperature. In the presence of 
Coulomb repulsion global phase separation becomes un- 
favorable and the system shows mesoscopic phase sepa- 
ration, where the size of charged regions is determined 
by the competition between ordering energy and the 
Coulomb energy. Using Monte-Carlo (MC) simulations 
we show that the system can form many different meso- 
scopic textures, such as clusters and stripes, depending 
only on the magnitude of the Coulomb repulsion com- 
pared to the anisotropic lattice attraction. Surprisingly, 
a feature arising from the anisotropy introduced by the 
Jahn- Teller interaction is that objects with even number 
of particles are found to be more stable than with odd 
number particles, which could be significant for super- 
conductivity when tunnelling is included|l8|. 

Let us consider the JT model Hamiltonian[l3T| . and 
take only the mode of B\ g symmetry: 
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where the Pauli matrix 03.1 describes the electronic dou- 
blet, g is a constant, and /(r) = (r 2 . — r 2 )/o(f) where 
fo(r) describes the effective range of the interaction^^. 

The model is reduced to a lattice gas model by using 
the adiabatic approximation for the phonon neld |l4Lll9| . 
The Hamiltonian in the pseudospin (5=1) representa- 
tion is given by: 

Hff-c = J2(~Vi(i -j)S?S? + V c (i - j)Q,Qj), (2) 

where Qi = {S[) 2 , V c {xa) — e 2 /eoam is 3D Coulomb 
potential, e is the charge of the electron, eo is the static 
dielectric constant and a is the effective lattice constant. 
S z = ±1 corresponds to the state with rii 2 = 1, rc.2,1 = 0> 
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and Sf = to ri\ = n 2 = 0. Simultaneous occupancy of 
both levels is excluded due to the large on-site Coulomb 
repulsion. The anisotropic short range attraction is then 
given by: 

V l (m)=g 2 /LoJ2f(i)f(m + i). (3) 

i 

A similar interaction can also be derived by considering 
the interaction of the electronic doublet with the strain 
of B\g symmetry, taking into account St.Venant's com- 
patibility conditions|lfil|. Anisotropic attraction caused 
by elasticity has the form[19j: 

^ (m) = -E exp( * k - m W 2 /^(k)) (4) 

here Aj are the components of the elastic modulus tensor, 
and U(k) = k i +s{ " Al / X 2)k i k i ■ Compared to (3), where 
the range of the interaction was defined by the coupling 
to optical phonons, the interaction (4) decays as 1/r 2 (in 
2D) at large distances. Since these attractive forces decay 
faster than the Coulomb repulsion at large distances, the 
net potential may have a minimum at short distances. 

Our goal is to study the model (2) at constant aver- 
age density of charged particles, n = jj-J^iQi- wnere 
N is the total number of sites. However, to clarify the 
physical picture we first consider a system with a fixed 
chemical potential by adding the term —fj, ^\ Qi to the 
Hamiltonian (2). 

Models such as (2), but in the absence of the long- 
range forces were previously studied on the basis of the 
molecular-field approximation |2flj . The mean- field equa- 
tions for particle density n and pseudospin magnetization 
M = i^i S? then have the form^: 

M = 2smh(2zV l M/k B T) 

exp(-u/k B T) + 2 cosh (2zViM/k B T) {1 

2 cosh (2zViM/k B T) 
exp(-u/k B T) + 2 cosh (2zViM/k B T) ( ' 

here z = 4 is the number of the nearest neighbours for a 
square lattice in 2D and k B is the Boltzman constant. A 
phase transition to an ordered state with finite M may 
be of either first or second order, depending on the value 
of /i. For the physically important case — 2zVi < fi < 
0, ordering occurs as a result of the first order phase 
transition. The two solutions of Eqs.(5,6) with M = 
and with M ^ correspond to two different minima of 
the free energy. The temperature of the phase transition 
T crit is determined by the condition: F(M = 0, fi, T) = 
F(M, fi, T) where M is the solution of Eq. (5) . When the 
number of particles is fixed (Eq.6), the system is unstable 
with respect to global phase separation below T cr i t . As 
a result, at fixed n two phases coexist with no = n(M = 
0,/i,T) and Um = n(M, /x, T), resulting in a liquid-gas- 
like phase diagram (Fig.l). 



To investigate the effects of the long range-forces, we 
performed MC simulations on the system (2). The sim- 
ulations were performed on a square lattice with dimen- 
sions up to L x L sites with 10 < L < 100 using a 
standard Metropolis al gori thmi c in combination with 
simulated annealing [22j |24j . At constant n one MC step 
included a single update for each site with nonzero Qi, 
where the trial move consisted from setting S z — at the 
site with nonzero Qi and S z — ±1 at a randomly selected 
site with zero Qi. A typical simulated annealing run 
consisted from a sequence of MC simulations at differ- 
ent temperatures. At each temperature the equilibration 
phase (10 3 — 10 6 MC steps) was followed by the averaging 
phase with the same or greater number of MC steps. Ob- 
servables were measured after each MC step during the 
averaging phase only. For L > 20 we observe virtually 
no dependence of the results on the system size. 

Comparing the MC results in absence of Coulomb re- 
pulsion shown by t cr it hi Fig. 1 with MF theory we find 
the usual reduction of t cr u due to fluctuations in 2D by 
a factor of ~ 2, . 

Next, we include the Coulomb interaction V c (r). We 
use open boundary conditions to avoid complications due 
to the long range Coulomb forces and ensure overall 
electroneutrality by adding a uniformly charged back- 
ground electrostatic potential (jellium) to Eq. (2). The 
short range potential vi(i) = VJ(i)eoa/e 2 was taken to be 
nonzero only for |i| < 2 and is therefore specified only 
for nearest, and next-nearest neighbours as vi(l,0) and 
vi(l, 1) respectively. 

The anisotropy of the short range potential has a pro- 
found influence on the particle ordering. We can see 
this if we fix vi (1,0) = —1, at a density n — 0.2 and 
vary the next-nearest neighbour potential U;(l, 1) in the 
range from —1 to 1. When uj(l, 1) < 0, the attraction 
is "ferrodistortive" in all directions, while for positive 
Vi(l, 1) > the interaction is " antiferrodistortive" along 
the diagonals. The resulting clustering and ordering of 
clusters at t = 0.04 is shown in Fig. 2a). As expected, a 
more symmetric attraction potential leads to the forma- 
tion of more symmetric clusters. On the other hand, for 
vi(l,l) = 1, the "antiferrodistortive" interaction along 
diagonals prevails, resulting in diagonal stripes. 

In the temperature region where clusters partially or- 
der the heat capacity (cl = d{E)i,/dT where E is the 
total energy) displays the peak at t co . The peak displays 
no scaling with L indicating that no long range ordering 
of clusters appears. Inspection of the particle distribu- 
tion snapshots at low temperatures (Fig. 2a) reveals that 
finite size domains form. Within domains the clusters are 
perfectly ordered. The domain wall dynamics seems to 
be much slower than our MC simulation timescale pre- 
venting domains to grow. The effective L is therefore 
limited by the domain size. This explains the absence of 
the scaling and clear evidence for a phase transition near 
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We now focus on the shape of the short range po- 
tential which promotes the formation of stripes shown 
in Fig. 2a). We set v t (l,0) = -1 and U;(l,l) = 
and study the density dependence. Since the inclu- 
sion of the Coulomb interaction completely suppresses 
the the first order phase transition at t cr jt, we mea- 
sure the nearest neighbor density correlation function 

9 P L = 4n(l-n)L" S|m|=l <Ei ~ Tt) (Qi - to 

detect clustering. Here () L represents the MC aver- 
age. We define a dimensionless temperature t c i — 
kBT c ieoa/e 2 as the characteristic crossover temperature 
related to the formation of clusters at which g p L rises to 
50% of its low temperature value. The dependence of t r j 
on the density n is shown in the phase diagram in Fig. 
1. Without Coulomb repulsion V c (r), t c i follows t crit , 
as expected. The addition of Coulomb repulsion V~c(r) 
results in a significant decrease of t c \ and suppression of 
clustering. At low densities we can estimate the onset for 
cluster formation by the temperature, to, at which g p L 
becomes positive. It is interesting to note that to almost 
coincides with the t cr n line at low n (Fig. 1). 

To illustrate this behaviour, in Fig. 2b) we show snap- 
shots of the calculated MC particle distributions at two 
different temperatures for different densities. The growth 
and ordering of clusters with decreasing temperature is 
clearly observed. At low n, the particles form mostly 
pairs with some short stripes. With further increasing 
density, quadruples gradually replace pairs, then long 
stripes appear, mixed with quadruples, etc.. At the high- 
est density, stripes prevail forming a labyrinth-like pat- 
tern. The density correlation function shows that the 
correlation length increases with doping, but long range 
order is never achieved (in contrast to the case without 
V c ). Note that while locally, there is no four-fold sym- 
metry, the overall correlation function still retains 4-fold 
symmetry. 

To get further insight in the cluster formation we mea- 
sured the cluster-size distribution. In Fig. 3 we show the 
temperature and density dependence of the cluster-size 
distribution function xl{]) — (N p (j)) L /{nL 2 ), where 
N p (j) is the total number of particles within clusters 
of size j. At the highest temperature xl{j) is close 
to the distribution expected for the random ordering. 
As the temperature is decreased, the number of larger 
clusters starts to increase at the expense of single parti- 
cles. Remarkably, as the temperature is further reduced, 
clusters of certain size start to prevail. This is clearly 
seen at higher densities (Fig. 3). Depending on the den- 
sity, the prevailing clusters are be pairs up to n « 0.2, 
quadruples for 0.1 < n < 0.3 etc.. We note that for a 
large range of w/(l,0), the system prefers clusters with 
an even number of particles. Odd particle-number clus- 
ters can also form, but have a much narrower parame- 
ter range of stability E^. The preference to certain clus- 
ter sizes becomes clearly apparent only at temperatures 
lower then t c ; , and the transition is not abrupt but grad- 



ual with the decreasing temperature. Similarly, with in- 
creasing density changes in textures also indicate a series 
of crossovers. 

The results of the MC simulation presented above al- 
low quite a general interpretation in terms of the kinetics 
of first order phase transitions |2^. Let us assume that 
a single cluster of order ed p hase with radius R appears. 
As was discussed in |l3j, UiJ |25j , the energy of the cluster 
is determined by three terms: e = —FirR 2 + awR + jR 3 . 
The first term is the energy gain due to the ordering 
phase transition where F is the energy difference between 
the two minima in the free energy density. The second 
term is the surface energy parameterized by a, and the 
third term is the Coulomb energy, parameterized by 7. 
If a < irF/Sj e has a well defined minimum at R — Ro 
corresponding to the optimal size of clusters in the sys- 
tem. Of course, these clusters are also interacting among 
themselves via Coulomb and strain forces, which leads 
to clusters ordering or freezing of cluster motion at low 
temperatures as shown by the MC simulations. 

We conclude that a model with only anisotropic JT 
strain and a long-range Coulomb interaction gives rise to 
a remarkably rich phase diagram including pairs, stripes 
and charge- and orbital- ordered phases, of clear rele- 
vance to functional oxides. The energy scale of the phe- 
nomena is defined by the parameters used in Hjt-c (2). 
For example, using the measured value eo — 40 |2q ] for 
La2Cu04, we estimate V c (l, 0) = 0.1 eV, which is also the 
typical energy scale of the "pseudogap" in the cuprates. 
The robust prevalence of the paired state in a wide re- 
gion of parameters (Fig. 3 c,d) is particularly interest- 
ing from the point of view of superconductivity. In con- 
trast to Bose condensation of mobile intersite bipolarons 
discussed by Alexandrov and Mott[13|, it has been sug- 
gested that pair tunnelling between objects such as shown 
in Fig. 2 can lead to an insulator-to-superconductor 
transition [18J. A similar situation occurs in mangan- 
ites and other oxides with the onset of a conductive state 
at the threshold of percolation, but different textures 
are expected to arise from the different magnitude (and 
anisotropy) of Vi(n), and static dielectric constant eq in 
the different materials [28). 
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Fig. 1. a) The phase diagram generated by Hjt 
(2) with, and without the Coulomb repulsion (CR). The 
dashed line is the MF critical temperature, while the full 
triangles (A) represent the MC critical temperature, t cr u, 
without CR. The open circles (o) represent t c i, without 
CR. The open triangles (A) represent t c i while the diag- 
onal crosses (x) represent the onset of clustering, t , in 
presence of CR. The cluster-ordering temperature (see 
text), t co , (also inch CR) is shown as crosses (+). The 
size of the symbols corresponds to the error bars, b) 
Typical temperature dependencies of the nearest neigh- 
bor density correlation function q l for n — 0.18 in ab- 
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Fig. 2. a) Snapshots of clusters ordering at t = 
0.04, n = 0.2 and w/(l,0) = -1 for different diago- 
nal vi(l, 1) (given in each figure). Grey and black dots 
represent particles clusters in state Sf — 1 and states 
= —1 respectively. The preference for even-particle- 
number clusters in certain cases is clearly observed, for 
example for 17/(1, 1) = —0.2. b) Snapshots of the particle 
distribution for two densities at two different tempera- 
tures t = 0.64 and t — 0.1 respectively. 




Fig. 3. The temperature dependence of the cluster- 
size distribution function xl{J) (for the smallest cluster 
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sizes) as a function of temperature at two different aver- and near t co (d) . The ranges of the density where pairs 
age densities n = 0.08 (a) and n = 0.18 (b). Xl(J) as a prevail are very clearly seen in (d). Error bars represent 
function of n at the temperature between to and t cl (c), the standard deviation. 



